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We adapted a multi-species escape model, developed for close-in extraso¬ 
lar planets, to calculate the escape rates of CH 4 and N 2 from Pluto. In the 
absence of escape, CH4 should overtake N 2 as the dominant species below 
the exobase. The CH4 prohle depends strongly on the escape rate, however, 
and the typical escape rates predicted for Pluto lead to a nearly constant mix¬ 
ing ratio of less than 1 % below the exobase. In this case the CH4 escape rate 
is only 5-10 % of the N 2 escape rate. Observations of the CH4 prohle by the 
New Horizons/ALICE spectrograph can constrain the CH4 escape rate and 
provide a nnique test for escape models. 


DRAFT 


July 20, 2017 


1:52cLm 


DRAFT 



KOSKINEN ET AL.: CH 4 ESCAPE FROM PLUTO 


X - 3 


1. Introduction 

Due to its low surface gravity, Pluto has an atmosphere that extends to several radii 
{Rp) and undergoes energy-limited escape. This means that a signihcant fraction of the 
solar EUV radiation absorbed by the upper atmosphere powers mass loss [Erwin et aL, 
2013; Zhu et aL, 2014], In this regard Pluto, despite its large distance from the Sun, 
is similar to many close-in extrasolar planets that undergo rapid, energy-limited escape 
[e.g., Yelle, 2004; Koskinen et aL, 2007, 2014], The same escape regime may also have 
shaped the composition of the early atmospheres of Venus, Earth and Mars [e.g., Zahnle 
and Easting, 1986; Hunten et aL, 1987; Tian et aL, 2008a, b], and is likely to be important 
on many potentially habitable rocky exoplanets to be studied in the future [e.g. Tarter 
et aL, 2007; Tian, 2009]. This makes the observations of Pluto’s upper atmosphere by 
the New Horizons/ALICE ultraviolet (UV) spectrograph especially interesting. 

The escape rate, that is constrained by the observed density and temperature prohles, 
is particularly important. It is the most pertinent test on the models, and a critical factor 
affecting the long-term evolution of the atmosphere. Most of the recent models, however, 
focus on N 2 , the dominant species. These models show that different escape mechanisms, 
hydrodynamic or Jeans escape, that are controlled by the upper boundary conditions, 
lead to different density and temperature prohles that nevertheless correspond roughly 
to the same mass loss rate [e.g., Erwin et aL, 2013]. This means that, depending on the 
model, the same observed density prohle can also be matched with different mass loss 
rates. Furthermore, the models indicate that the density prohles in the upper atmosphere 
are sensitive to thermal structure and radiative balance in the lower atmosphere that are 
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not directly related to the escape rate [Zhu et ai, 2014], These factors complicate the use 
of the N 2 prohle to constrain the escape rate. 

The ALICE spectrograph covers the wavelength range from 50 nm to 190 nm [Stern 
et ai, 2005], including the ionization and dissociation continua of N 2 and CH4, as well 
as the N 2 electronic bands, that probe the upper atmosphere above r Ri 1.5 Rp. Thus 
the solar occultation measurements that will be obtained during the encounter provide 
simultaneous coverage of N 2 and CH4 at radial distances of about 1.5-3 Rp. The density 
prohle of CH4 is controlled by both diffusion and the escape hux, and it can therefore 
provide further constraints on the escape rate. In this work we use a multi-species escape 
model, developed originally for extrasolar planets [Koskinen et al, 2013a, b], to calculate 
the escape rates of N 2 and CH4 from Pluto at the conditions expected during the New 
Horizons encounter. 

We note that early models of escape from Pluto assumed that CH4 is the dominant 
species in the upper atmosphere, based on the argument that diffusive separation should 
lead CH4 to overtake N 2 [McNutt, 1989; Hubbard et al, 1990]. In contrast, later work by 
Krasnopolsky and Cruikshank [1999], that included the effect of escape on CH4 and several 
other minor species in addition to diffusion, indicated that N 2 is the dominant species at all 
altitudes. However, they also showed that the abundance of CH4 is potentially signihcant 
near the exobase. The purpose of this work is to revise the predictions for CH4 escape 
in light of the recent progress in models of atmospheric escape for Pluto [Tucker et al, 
2012; Erwin et al, 2013; Zhu et al, 2014], and to study the sensitivity of the CH4 density 
prohle to the escape rate, diffusion and surface mixing ratio. 
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2. Methods 

A detailed description of our escape model is given in Koskinen et al. [2013a, b]. The 
model solves the time-dependent equations of continuity, momentum and energy in the 
vertical direction. We note that a time-dependent model passes smoothly through the 
critical point in transonic escape, although this is not important here because escape 
from Pluto is subsonic. We solve common momentum and energy equations for the 
bulk flow, together with separate continuity equations for the different species i.e., N 2 , 

CH4 and CO. The velocities of the these species are solved self-consistently by using the 
diffusion approximation [e.g., Garcia Munoz, 2007]. The model is capable of including 
photochemistry in the upper atmosphere but this option is not used here because its 
influence on the density prohles should be negligible. 

In order to calculate the solar EUV heating rate in the upper atmosphere, we used a 
Level 3 spectrum of the Sun from the TIMED/SEE archive (http;//lasp.colorado.edu/lisird/see/) 
generated for 28 April, 2015 at a resolution of 1 nm. We calculated the absorption of ra¬ 
diation self-consistently in the model, and assumed EUV heating efficiencies of 25 % and 
50 % for N 2 and CH4, respectively [Krasnopolsky, 1999]. This calculation differs from 
Zhu et al. [2014] in that we used the full solar spectrum to calculate the EUV heating 
rate instead of only 3 broad wavelength bins. We also used the parameterization given 
by Strobel [2008] to include infrared (IR) cooling from rotational lines of CO. We did not 
include heating and cooling by the near-IR bands of CH4 that should be negligible in the 
upper atmosphere above our lower boundary at po = 10“^ bar. 
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Boundary conditions are of critical importance to the results. At the lower boundary 
we set the temperature to Tq = 103.7 K and altitude to Zq = 340.5 km, in agreement 
with the reference model of Zhu et al. [2014]. For the composition, we varied the volume 
mixing ratio of CH 4 by adopting surface values of 0.3 %, 0.44 % and 0.6 % for different 
simulations, based on the range derived from VLT/CRIRES observations [Lellouch et a/., 
2015]. In line with Zhu et al. [2014], we assumed a surface mixing ratio of 5 x 10“^ for 
CO in all of our simulations. The eddy mixing rate is currently unknown, and we adopted 
values of 10, 100, and 1000 m^ s“^ for the eddy diffusion coefficient K^z- The molecular 
diffusion coefficients were taken from the compilation of Marrero and Mason [1972]. Our 
reference model assumes a CH4 surface mixing ratio of 0.44 % and Kzz = 100 m^ s“^. 
The CH4-N2 diffusion coefficient in the reference model varies from 34 m^ s“^ at the lower 
boundary to 8.5 x 10^ m^ s“^ at the exobase. 

We can use the above information to verify that the CH4 profile is not significantly 
affected by photochemistry. The photolysis of CH4 is dominated by the solar Lyman a 
photon flux, which is 3.8 x 10^^ m“^ s“^ at 32.9 AU based on the TIMED/SEE spectrum. 
The CH4 cross section at Lyman a is 1.85 x 10“^^ m^ and the subsolar optical depth 
based on our model reaches unity near 1.4 Rp. The subsolar (global minimum) photolysis 
timescale at 1.4 Rp is therefore 3.9 x 10® s. The CH 4 -N 2 diffusion coefficient at that level 
is Dst = 1.7 X 10^ m^ s“^ and the scale height is 92 km. Thus the subsolar photolysis 
timescale is longer than the diffusion timescale of = H'^/Dgt = 5 x 10^ s. This difference 
increases rapidly with altitude and near the effective heating peak of about 1.8 Rp (see 
Section 3.3) the photolysis timescale is 50 times longer than the diffusion timescale. 
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At the upper boundary, we used the modihed Jeans (Type-I) boundary conditions with 
an enhancement factor of 1.5 [Zhu et ai, 2014], The npper bonndary is placed at the 
exobase (where the Knudsen number Kn = 1), the altitude of which is automatically 
adjusted by our model while the simulations approach steady state. We used the same 
boundary condition for the temperature gradient as Zhu et al. [2014]. Their model, how¬ 
ever, did not consider CO and CH 4 separately whereas onr model reqnires additional 
bonndary conditions for the velocities and partial pressnre gradients of the individnal 
species. 

We set the velocities of the different species at the exobase to the Jeans effusion velocity: 


Wsj = Ti 


I kT 
2nm. 


(1 -h AJexp {-\s) 


( 1 ) 


where is the escape parameter for a species with mass rUg and T = 1.5. The pressnre 
gradient for species s at the upper boundary can then be derived from the diffusion 
eqnation [e.g., Schunk and Nagy, 2000], and it is given by: 


_ ^1^ _ V —-lV — 

Ps dr mpdr Dgt 


( 2 ) 


t^S t^S 

where p = J2sPs, ur is the mean molecular weight, Xt is the volume mixing ratio and Dgt 
is the mutnal diffusion coefficient for species s and t. We note that eddy diffusion does 
not need to be included in the npper bonndary conditions because K^z « Dg at the 
exobase. 


3. Results 
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3.1. CH 4 mixing ratios 

The CH4 abundance profile in Pluto’s atmosphere is controlled primarily by escape, 
with limited sensitivity to our range of surface mixing ratios or K^z- This is illustrated by 
Figure 1 that shows the CH 4 prohles for 5 different simulations. The solid lines show results 
for models with K^z = 10-1000 m^ s“^ and a surface mixing ratio of 0.44 %, including 
the reference model, while the dashed and dashed-triple-dotted lines show models with 
Kzz = 100 m^ s“^ and CH4 surface mixing ratios of 0.3 % and 0.6 %, respectively. The 
global escape rate in these models, which increases with the surface mixing ratio of CH4, is 
1.8-2.3 X 10^’^ s“^. This is the sum of the escape rates of the individual species (hereafter, 
the total escape rate). For comparison, the dotted and dash-dotted lines show the zero 
escape and reduced escape CH4 prohles, respectively, based on the T-P prohle in the 
reference model. The reduced total escape rate here is 2.9 x 10^® s“^, based on a reduced 
Jeans enhancement factor of P = 0.15. 

In the absence of escape CH4 would become the dominant species in the reference model 
around r = 4.4 Rp i.e., well below the exobase, which is located at r = 5-6 Rp. With 
escape included, however, the CH4 mixing ratio is less than 1 % at all altitudes. Near 
r = 3 Rp where solar occupations in the EUV become sensitive to absorption by CH4 (see 
Section 3.5), the mixing ratio in the reference model is only 0.46 % whereas in the zero 
escape case it is almost 50 times larger at 21.4 %. An escape rate of about 2 x 10^^ s“^ 
therefore leads to a nearly constant mixing ratio of CH4 below r = 3 Rp. 

The strong dependence of the CH4 mixing ratios on escape means that observations of 
the CH4 prohle provide a means to constrain the escape rate and thus the hrst opportunity 
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to partially validate a hydrodynamic escape model for a planetary atmosphere. The escape 
rate predicted by the current models is near the saturation value i.e., the value that leads 
to a nearly constant mixing ratio, which is close to the surface mixing ratio, below 3 Rp. 
As shown by Figure 1, a lower total escape rate of 2.9 x 10^® s“^ with the reference model 
T-P prohle leads to a higher CH 4 mixing ratio of 1.5 % at r = 3 Rp. We note that the 
CH4 prohles in Figure 1 are practically independent of Kzz-, which further enhances their 
potential to constrain the escape rate. 

We note that our results differ from Krasnopolsky and Cruikshank [1999] who also 
studied the diffusion and escape of CH4 on Pluto. Their model, which included an N 2 
escape rate of 2.6 x 10^^ s“^, predicted that diffusive separation of CH4 takes place near 
r = 2300 km, and that the mixing ratio of CH4 is about 19 % at r = 3500 km. In our 
models the mixing ratio of CH4 at r = 3500 km varies between 0.32 % and 0.61 %. The 
differences between our model and the earlier work are difficult to identify exactly. The 
exobase, however, was typically at a lower altitude in the model of Krasnopolsky and 
Cruikshank [1999], and their temperatures were generally cooler. Our temperatures are 
signihcantly warmer due to the upper boundary condition based on the Jeans energy flux. 

3.2. Temperature profiles 

The temperature prohles based on our models are shown by Figure 2. The differences 
in the predicted temperatures are within 5 K. The peak temperature increases with the 
abundance of CH4, while the temperature near the exobase decreases slightly as the peak 
temperature increases. This is caused by more efficient absorption of Lyman a radiation 
by CH4, which leads to both higher peak temperatures and a slightly higher escape rate. 


DRAFT 


July 20, 2017, l:52cLm 


DRAFT 



X - 10 


KOSKINEN ET AL.: CH 4 ESCAPE FROM PLUTO 


As a result, adiabatic cooling near the exobase is also more efficient, leading to a slightly 
lower temperature near the exobase. Due to the nearly constant mixing ratio of less than 
1 % for CH4 in the model, however, the inclusion of diffusion does not signihcantly affect 
the results of Zhu et al. [2014] who assumed that CH 4 is uniformly mixed with the surface 
value of 0.44 %. Our temperature prohles do differ slightly from their predictions, but 
this is mostly due to the CO cooling parameterization in our work that leads to lower 
peak temperatures. 

3.3. Energy balance and CH4 escape 

Our reference model predicts a total escape rate of 2.1 x 10^^ s“^. This value is es¬ 
sentially identical to the previous estimates of the escape rate from Pluto’s atmosphere. 
We hnd that N 2 and CO escape from the exobase with an outflow velocity of 1.75 m s“^, 
while the CH4 outflow velocity is much higher, about 18 m s“^. Due to the relatively 
low abundance of CH4, however, the effective bulk outflow velocity i.e., w = J2s PsWg/p, 
is still only 1.82 m s“^. The predicted escape rate of CH4 in our models is 1-2 x 10^® 
s“\ or about 5-10 % of the N 2 escape rate. These values are an order of magnitude lower 
than the diffusion limit of 2 x 10^^ s“^ [Strobel, 2008]. Thus our reference model agrees 
with recent work in that escape from Pluto is strongly subsonic, and that the escape of 
N 2 dominates [e.g., Erwin et al, 2013; Zhu et al, 2014]. 

These models also show that mass loss from Pluto is energy-limited. This statement 
needs to be qualihed, because of the efficiency factors that affect this conclusion. The 
energy limited mass loss rate based on globally averaged solar flux is: 

M = (3) 
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where Fe is the integrated solar flux, in this case at wavelengths of 1-145 nm, rjE is the 
mass loss efficiency, is the heating efficiency, te is the effective radius of the heating 
peak and A<l> is the gravitational potential difference that the escaping particles have to 
overcome. 

At the distance of 32.9 AU, the integrated solar flux is Fe = 1.18 x 10“^ W m“^, based 
on the TIMED/SEE spectrum. The global, effective heat flux in our reference model is 
oieFe = 3.01 X 10“® W m“^, which is based on the net heating rate that includes CO 
cooling, and te = 1.78 Rp. Thus the heating efficiency is aE ~ 0.26, and the energy- 
limited mass loss rate based on this efficiency would be 103 kg s“^. The mass loss rate 
predicted by the reference model, on the other hand, is 96 kg s“^, implying that the mass 
loss efficiency is tje ~ 0.93. This means that, at least according to the models, escape 
from Pluto is energy-limited i.e., tje —)■ 1, in analogy to many close-in extrasolar giant 
planets [e.g., Koskinen et ai, 2014], and in contrast to the present planetary atmospheres 
in the solar system. 

3.4. Boundary conditions 

Boundary conditions for escape models have been discussed extensively in recent liter¬ 
ature on Pluto, but so far this discussion has received only cursory attention elsewhere. 
For example, many of the existing time-dependent exoplanet models rely on the so-called 
outflow boundary conditions i.e., extrapolation of the density, temperature and velocity 
with a constant slope at the upper boundary [e.g., Tian et ai, 2005; Garcia Munoz, 2007]. 
These boundary conditions lead to unphysical outcomes on Pluto, which has an exobase 
well below the altitude of the implied sonic point. In theory, the outflow boundary condi- 
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tions may be valid above the sonic point, although this assertion probably deserves further 
attention in future studies. 

Other exoplanet models have imposed Jeans or modihed Jeans outflow velocities at the 
upper boundary with extrapolated temperatures, either with a zero gradient or a constant 
slope [e.g., Tian et a/., 2008a, b; Koskinen et al, 2014], This approach appears to be valid 
only if the atmosphere is rendered isothermal by conduction at the exobase, in which 
case the constant slope also reduces to the isothermal boundary condition and produces 
acceptable results. On Pluto, extrapolating the temperature with a constant slope leads 
to an almost adiabatic decrease in temperature with altitude, and a much cooler exobase. 
As a result, the conductive heat flux becomes larger than the Jeans energy flux. Zero 
gradient, on the other hand, leads to an overestimated temperature near the exobase. 
In general, the validity of the boundary conditions in the present models, and thus the 
escape solution, can be tested by matching the models with both the N 2 and CH4 density 
prohles that are, at least in theory, retrievable from the EUV solar occultations. 

3.5. Observations 

The model calculations show that the CH4 profile can be used as an indicator of the 
escape rate; however, it is important to have coincident measurements of N 2 because this 
enables calculation of the CH4 mixing ratio and the atmospheric temperature. We note 
that stellar and solar occultations observed by the Cassini/UVIS instrument in its EUV 
channel have been used to retrieve coincident density prohles of N 2 and CH4 in Titan’s 
upper atmosphere for the same range of slant column densities that is relevant for escape 
models on Pluto. In particular, Capalbo et al. [2013] used the 58.4 nm and 63 nm solar 


DRAFT 


July 20, 2017, l:52cLm 


DRAFT 



KOSKINEN ET AL.: CH 4 ESCAPE FROM PLUTO 


X - 13 


emission lines in the ionization continuum of N 2 , together with a wavelength bin around 
108.5 nm in the ionization continuum of CH 4 , to retrieve these density prohles from a 
solar occultation. In addition, Kammer et al. [2013] used EUV stellar occultations at 
91.1-110 nm, probing the N 2 electronic band system, for the same purpose. 

For example, we used our models to calculate predicted transmission for Pluto in the 
108.5 nm bin, which includes a group of solar N II and He II emission lines [Curdt et al, 
2001]. In this region absorption is by CH4 and the cross section is roughly constant 
with wavelength [Kameta et al, 2002], allowing for a simple transmission calculation 
instead of a more complex forward model. The results are shown in Figure 3, which 
indicates that the zero escape case is clearly distinguishable from our reference model. 
The difference in transmission between the reference model and the reduced escape model 
is less pronounced, peaking at about 0.09 around 1.9 Rp where transmission is about 0.7. 
This difference is detectable with a signal-to-noise ratio (S/N) of 8 or higher. 

We also used our reference model to calculate transmission in the 63 nm wavelength 
bin, which includes the solar O V line that is absorbed predominantly by N 2 , with a small 
contribution from CH 4 . The resulting light curve is limited to altitudes above 2.5 Rp 
(Figure 3). At lower altitudes the N 2 electronic band system (80-100 nm) can be used to 
retrieve the N 2 prohle. To illustrate this, we calculated transmission at 89-90.5 nm in the 
solar Lyman continuum, in order to avoid uncertainties associated with the width of the 
solar emission lines. Some of these lines are absorbed by the N 2 bands while absorption 
by CH 4 features strongly between the bands. A high resolution forward model is therefore 
required to calculate transmission, even when instrument broadening is signihcant. In our 
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model we used a SUMER spectrum of the Sun [Curdt et ai, 2001] and a high resolution 
N 2 band cross section [Lewis et ai, 2008]. In order to include the point spread function 
(PSF) for ALICE, we used a Gaussian with a FWHM of 0.4 nm [Stern et ai, 2005]. The 
resulting occultation light curves in Figure 3 show that CH4 can be measured from 1.6 Rp 
to 2-4 Rp and N 2 can be measured from 1.8 Rp to 4 Rp, assuming that densities can be 
retrieved for transmissions between 0.1 and 0.9. The precise altitude range will, of course, 
depend on the noise characteristics of the data set. 

4. Conclusions 

We adapted a multi-species model of hydrodynamic escape to calculate the escape rates 
of N2 and CH4 from Pluto, and to quantify the effect of escape on the CH4 prohle. The CH4 
prohle plays a critical role in the heating and photochemistry of the upper atmosphere, 
and under diffusive separation it should overtake N 2 as the dominant species below the 
exobase. As we showed, it can also be used to constrain the escape rate. This provides a 
unique opportunity to test models of hydrodynamic escape in an actual atmosphere. 

The total escape rate predicted by our reference model is 2.1 x 10^^ s“^, in agreement 
with previous models. In terms of energy-limited escape, this corresponds to ~24 % of 
the incident solar UV flux at 1-145 nm with an effective heating peak around 1.78 Rp. 
The CH 4 escape rate in the model is about 5-10 % of the N 2 escape rate. We hnd that 
the mixing ratio of CH 4 in the upper atmosphere depends strongly on the escape rate and 
is almost independent of Kzz- The total escape rate of 2.1 x 10^'^ s“^ leads to a nearly 
constant mixing ratio that is close to the surface value below r = 3 Rp whereas lower escape 
rates lead to substantially higher mixing ratios. Given an independent measurement of 
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the surface mixing ratio, the CH 4 and N 2 prohles at r = 1.8-3 Rp that are retrievable 
from EUV occultations can therefore constrain the escape rate. 
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CH4 mixing ratio 


Figure 1. The mixing ratio of CH 4 in different simulations. The solid lines (that are indistin¬ 
guishable) are based on a surface mixing ratio of 0.44 % and K^z = 10-1000 m^ s“^. The dashed 
and dash-triple-dotted lines are based on surface mixing ratios of 0.3 % and 0.6 %, respectively, 
with Kzz = 100 m^ s“^. The dotted and dash-dotted lines show the zero escape and reduced 
escape rate prohles, respectively, based on our reference model temperature-pressure prohle, a 
surface mixing ratio of 0.44 % and K^z = 100 m^ s“^. In each case, the upper boundary is the 
exobase. 
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Figure 2. Temperature profiles based on simulations with different surface mixing ratios 
and Kzzi corresponding to the CH 4 mixing ratio prohles in Figure 1. The solid lines (that are 
indistinguishable) show results for models with a CH4 surface mixing ratio 0.44 % and K^z = 10- 
1000 m^ s“T The dashed lines and dash-triple-dotted lines show results based on CH4 surface 
mixing ratios of 0.3 % and 0.6 %, respectively, with Kzz = 100 m^ s“T 
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Figure 3. Predicted transmission for different EUV wavelength bins. The red lines show 
transmission in the 108.5 (108-109) nm bin based on the CH 4 profiles in the reference (solid 
line), zero escape (dotted line) and rednced escape (dash-dotted line) models (see Fignre 1). The 
bine lines show transmission in the 63 (62.5-63.5) nm bin (dashed line) and 90 (89-90.5) nm bin 
(solid line) based on the reference model N 2 and CH 4 profiles. 
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